Langevin Dynamics of the vortex matter two-stage melting transition in 
Bi 2 Sr2CaCu 2 8+( 5 in the presence of straight and of tilted columnar defects 
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In this paper we use London Langevin molecular dynamics simulations to investigate the vor- 
tex matter melting transition in the highly anisotropic high-temperature superconductor material 
Bi2Sr2CaCu208+i in the presence of low concentration of columnar defects (CDs). We reproduce 
with further details our previous results obtained by using Multilevel Monte Carlo simulations that 
showed that the melting of the nanocrystalline vortex matter occurs in two stages: a first stage melt- 
ing into nanoliquid vortex matter and a second stage derealization transition into a homogeneous 
liquid. Furthermore, we report on new dynamical measurements in the presence of a current that 
identifies clearly the irreversibility line and the second stage derealization transition. In addition 
to CDs aligned along the c-axis we also simulate the case of tilted CDs which are aligned at an angle 
with respect to the applied magnetic field. Results for CDs tilted by 45° with respect to c-axis show 
that the locations of the melting and derealization transitions are not affected by the tilt when 
the ratio of flux lines to CDs remains constant. On the other hand we argue that some dynamical 
properties and in particular the position of the irreversibility line should be affected. 

PACS numbers: 74.25.Qt, 74.25.Ha, 74.25. Dw, 74.25.Bt 



I. INTRODUCTION 

Since their discovery during the 1980s, much progress 
has been made on the study of high-temperature super- 
conductors. These materials are classified as type II su- 
perconductors, for which there is only a partial Meissner 
effect when the external magnetic field is in the range 
H c i < H < H c -2r*2£. Inside the superconductors, the 
magnetic field survives in the form of quantized flux lines 
(FLs), called vortices, each with a quantum unit of flux 
4>q = hc/2e. The vortices form a periodic hexagonal lat- 
tice (Abrikosov lattice) at low temperatures and melt into 
a vortex liquid at higher temperatures^&LIui. 

High-temperature superconductors are ceramic mate- 
rials which have the common structure of layered super- 
conducting copper-oxide planes. They are anisotropic 
materials, with anisotropy defined as 7 = y / m z /m±, 
where m z and m± are the effective masses of electrons 
moving perpendicular to the superconducting planes and 
along to them respectively. One of the commonly studied 
materials Bi 2 Si'2CaCu20g + 5, known as BSCCO, is highly 
anisotropic. It is estimated to has an anisotropy in the 
range of 300-500 for optimally doped samples. 

In high anisotropic materials, vortices can be described 
as stacks of pancake vortice o 10 ! 11 . The "pancakes" are 
centered at their corresponding planes and can only move 
within their own planes. The interaction between the 
pancakes consists of two parts: the electromagnetic in- 
teraction and the Josephson interaction. 

The electromagnetic interaction originates from the 
screening currents that arise in the same plane where 
a pancake resides as well as in more distant planes. It is 
repulsive among pancakes in the same plane and attrac- 
tive for those in different plane o 10 ' 12 . This interaction 
exists even in the case where the superconducting planes 
are completely decoupled, so no current can flow along 



the c axis (perpendicular to the planes) of the sample . 

The Josephson interactio n 2 i 12 i 13 originates from the 
Josephson current flowing between the superconduct- 
ing planes, which is proportional to the sine of the 
gauge invariant phase difference between two planes 
A(p + (2ir / (f>o) J ds ■ A, where ip is the phase of the super- 
conducting order parameter and A is the vector poten- 
tial. When two pancakes belonging to the same stack and 
residing in adjacent planes move away from each other, 
the phase difference that originates causes a Josephson 
current to begin flowing between the planes. This re- 
sults in an attractive interaction between pancakes that 
for distances small compared to A = jd is approximately 
quadratic* 2 -^ in the distance. Here we denoted by d the 
inter-plane separation. When the two adjacent pancakes 
are separated by a distance larger than A, a "Joseph- 
son string" is formed, whose energy is proportional to its 
length 1 ^! 4 .. 

Recently the effect of columnar defects (CDs) 
on the melting transition was investigated both 
exp eriment all y 1 5 1 1 6 1 1 7 and numericall y 18 ^ 9 ' 20 ! 21 . Colum- 
nar defects can be induced artificially by irradiating the 
sample with high energy heavy ions, like I GeV lead 
ions. The heat released by the ions create tubes of 
damaged, non-conducting material. It is customary to 
measure the density of CDs by the value of the "match- 
ing field" by = n c d<po, where n ct t is the number of 
CDs per unit area. In the limit B <C the vortex 
solid becomes a Bose glass as most FLs are trapped by 
CDs located at random positions 2 ^ 22 ! 23 ! 24 . On the other 
hand our interest is in the case B^ < B, in which there 
are more FLs than CDs. In this case the picture that 
emerged from experiment a 15 ' 16 ! 17 and from numerical 
investigation a 18 ' 19 ! 20 ' 21 and analytical calculation o 21 ' 25 is 
that of the "crystallites in the pores" , also referred to as 
the "porous" vortex solid. At low temperatures a skele- 
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ton (or matrix) of vortices localized on CDs is formed, 
while the excess (interstitial) FLs try to form hexagonal 
crystallites in the lacuna between CDs. Thus this phase 
has a short ranged translational order, which extends to 
a distance of the order of a typical pore size. 

As the temperature is increased, and if the magnetic 
field is large enough, this heterogeneous structure melts 
in two stages: first the crystallites in the pores melt into 
a nanoliquid while the skeleton remain intact, and subse- 
quently the skeleton melts and the liquid becomes homo- 
geneous. When the magnetic field is lowered these two 
transitions coincide into one. This is usually associated 
with a kink in the first order melting line. The first stage 
melting of the crystallites is observed in the experiments 
as a step in the equilibrium local magnetizatio n 26 ' 27 i 28 
and in the simulations as a sharp increase in the trans- 
verse fluctuations of the FLs, among other signatures. 
The second transition was not observed in the exper- 
iments as a similar jump in the magnetization or an- 
other equilibrium property. It was observed in transport 
measurement^ involving transport currents with alter- 
nating polarity. This establishes the transition as dynam- 
ical in nature. It remains to be established that this is 
also an equilibrium thermodynamic transition. In order 
to support this assertion it was arguecUi that the second 
"derealization" transition is associated with the restora- 
tion of a broken longitudinal gauge symmetry. However, 
at this time it cannot be ruled out that the second tran- 
sition is only a crossover associated with a gradual equi- 
librium change over a finite range of temperatures. 

Recently one of us published a short paper— report- 
ing on multilevel Monte Carlo simulations of the vortex 
matter in the highly anisotropic high-temperature super- 
conductor BSCCO. We introduced low concentration of 
columnar defects satisfying B^ < B. Both the electro- 
magnetic and Josephson interactions among pancake vor- 
tices were included. The nanocrystalline, nanoliquid, and 
homogeneous liquid phases were identified in agreement 
with experiments. We observed the two-step melting pro- 
cess and also noted an enhancement of the structure fac- 
tor just prior to the melting transition. We also proposed 
a simple theoretical model to explain some of the find- 
ings. In this paper we would like to give more details 
of the simulation results but we use a different method 
than in the previous work, i.e. we use molecular dynam- 
ics instead of multilevel Monte carlo. This enables us to 
carry out dynamical measurements that clearly identifies 
the second melting transition (derealization) . We also 
report in this paper simulations of tilted columnar de- 
fects. A short summary of our findings for tilted defects 
has recently been submitted 2 ^. 

Molecular dynamics (MD) is a powerful tool for sim- 
ulations of physical systems and it often serves as an 
alternative to Monte Carlo (MC) simulations. Its ad- 
vantage is that it can be used to investigate the real 
dynamics of the system as opposed to MC simulations 
that are used for obtaining equilibrium properties. How- 
ever MD simulations could be plagued by the absence of 



ergodicity when applied to systems represented by path 
integrals 3 -^ and there is also the problem of implement- 
ing permutations for the case of identical particles like 
Bosons. The problem of ergodicity is really not much 
of an issue for Langevin simulations since the thermal 
noise helps the system explore the configuration space 
and it can be shown by using the corresponding Fokker- 
Planck equation that equilibrium is reached in the long 
time limit. For flux-lines we have found a way to im- 
plement "permutations" in the MD simulations by flux 
cutting and recombining and this was explained in a re- 
cent publication 31 which reported on MD simulations of 
the melting transition in pure BSCCO without colum- 
nar or point defects. We were also able to implement 
periodic boundary conditions in all directions (including 
the z direction), and to include the in and out of plane 
electromagnetic interaction as well as the Josephson in- 
teraction using the new approximation we have recently 
obtained^ 2 -. Results for the first-order melting transitions 
in BSCCO for fields of 100-200 gauss were found to be in 
an excellent agreement with the results of our previous 
multilevel Monte Carlo simulations 3 - 3 - including the pro- 
liferation of non-simple loops corresponding to flux en- 
tanglement above the melting transition. In this paper 
we extend the MD simulations to systems with columnar 
defects, both for the case that the columns are paral- 
lel to the c-axis and for the case when the columns are 
tilted at an angle 9 with respect to the c-axis. The ex- 
ternal field is taken to be aligned along the c-axis. We 
would like to check if the location of the melting tran- 
sition in the B-T plane is sensitive to the tilting angle 
of the columnar defects. Recently experiments were car- 
ried out to investigate this problem 2 ^. The conclusions 
were that while static thermodynamical signatures like 
the location of the melting and derealization lines do 
not depend (within the experimental errors) on the an- 
gle between the applied field and the columnar defects, 
non-equilibrium properties like the position of the irre- 
versibility line do depend on the angle. 

II. SIMULATION METHOD 
A. method 

In this paper we use the London-Langevin molecular 
dynamics simulation method. The details of the simu- 
lation method, for a system without columnar defects is 
given in our previous paper—. Here we summarize the 
main ingredients of the method. The equation of motion 
for the m'th pancake vortex is 

d V^p = - V m V({r n }) + f L + Cn(t). (2.1) 

The pancake label m stands actually for two indices (i,p) 
where p is the plane label and i is the pancake label in 
that plane. The position R m is a two component vec- 
tor in the plane. We use the over-damped model for 
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vortex motion in which the velocity of the vortex is pro- 
portional to the applied force and 77 is the viscous drag 
coefficient per unit length given by the Bardeen-Stephen 
expression^!. In Eq. (|2.I[) d is the interlayer spacing be- 
tween CuC>2 planes that is taken to be equal to the width 
of the pancake vortex. V is the potential energy depend- 
ing on the position of all pancakes and includes both the 
magnetic energy and Josephson energy, and the interac- 
tion between a pancake and the columnar defects. The 
force is minus the gradient of the potential energy with 
respect to the position of the m'th pancake, fa is a 
driving force (if present), for example the Lorentz force 
induced by a current. £ m is a white thermal noise term 
which satisfies 

(G(*)Cf (*')> = 2kT V d S a0 S mn S(t - t'). (2.2) 



In Eq. (|2.2p a and (i refer to the x and y components 
of the vector £ and m and n are pancake labels, k is 
Boltzmann's constant. In the simulations we use a dis- 
crete form of the Dirac delta-function. We measure dis- 



tances in units of ao = y 2<f>o/B-\/3 where B is the mag- 
netic field. We measure energy in units of eod where 
eo(T) = (4>o/4:TtX) 2 is the basic energy scale per unit 
length and A is the penetration depth. We measure 
time in units of rya§/eo. For example for T = 60K 
and B = 100G, a « 4887 A, e d « 4.685 x 10~ 14 erg 
1=3 339.5 K/k. Ref. j3f| quotes a value for 77 for a single 
crystal BSCCO of around 1 x 10~ 7 g/(cm s). Based on 
this value the time unit is about 0.765 ns. These values 
change when the temperature and field are varied. The 
simulation cell was chosen to have a rectangular cross 
section of size L\ x L2 — ao y/Nji x ao\/3Nfi/2 where 
Nfi is the number of flux lines (number of pancake vor- 
tices in each plane). We usually worked with 36 flux lines 
but to check for finite size effect we also used 64 FLs in 
some simulations. The aspect ratio of the cell was chosen 
to accommodate a triangular lattice without distortion, 
such that each triangle is equilateral. In the z-direction 
we take N p layers of width d each, where in practice we 
have chosen N p to be 36 or 200, as will be indicated be- 
low. Thus the number of pancakes used is from 1296 up 
to 7200. Note that our system is effectively bigger since 
every pancake interacts not only with all other pancakes 
in the simulation cell but with infinitely many images of 
these pancakes in neighboring cells in all directions. 



B. 



interactions 



The electromagnetic (EM) interaction occurs between 
any pair of pancakes. To produce more realistic results 
for the small system simulated we have implemented peri- 
odic boundary conditions (PBC) in all directions includ- 
ing the z-direction. Periodic boundary conditions mean 
that every pancake interact not only with the actual pan- 
cakes in the simulation cell but also with all their images 
in other cells which are part of an infinite periodic array. 



Each image of a pancake is located at the same position 
in the corresponding cell as the original pancake in the 
simulation cell. The details are explained in Ref. [3l| . 
The final expression we used for the pair energy of two 
pancakes with fully implemented PBC is 
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for pancakes in different planes separated by a vertical 
distance dAp and a lateral separation R = (x.y), and 
similarly 
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for pancakes in the same plane. The functions 
Go(R/X, L\/X) and Gqc{x/ Li, yjL-i) represent the peri- 
odic screened and bare Greens functions for the Coulomb 
interactions in 2D and are given in the Appendix of 
Ref. [3l[ . The factor f m comes from the summation of the 
EM interaction between a pancake and another pancake 
in the same cell and all its images along the z-direction. 
It is given by 



/ ro (Ap)= J2 exp(-|Ap + iVj,%) 



2=-c 



exp(-|Ap| M ) + exp((|Ap|-Ag M ) 
1 - exp(-N p fi) 



(2.5) 



where we put /_t = d/X. The dependence of f m on Ap 
is rather weak. The approximations are valid provided 
dN p < X. 

We have also implemented the Josephson interaction 
among adjacent pancakes belonging to the same 
The approximate formula we use is based on our recent 
paper that derives a numerical solution to the nonlinear 
sine Gordon equation in two dimensions^. The formula 
interpolates between an eod(R/A) 2 ln(A/i?) dependence 
for small R to e Q dR/A dependence for large R, where R 
is the lateral separation of the pancakes, A = jd, 7 is 
the anisotropy and d is the inter-plane separation. This 
interpolation is an modification of a similar interpolation 
formula previously used by Ryu, Doniach, Deutscher and 
Kapitulnik 39 . The Josephson interaction diminishes with 
increasing anisotropy, and we choose the anisotropy pa- 
rameter to be 375 to obtain reasonable agreement with 
recent experiment on pristine melting. In this work we 
neglect Josephson interaction between pancakes belong- 
ing to different stacks because on the average those pan- 
cakes are farther away. In addition we view the stack 
as the collection of pancakes carrying the unit quantum 
of flux from one side of the sample to the other. Since 
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the josephson interaction is proportional to the gauge in- 
variant phase difference Aip + (27r/0 o ) J ds ■ A, the sec- 
ond term involving the vector potential will be larger for 
pancakes belonging to the same stack than to pancakes 
belonging to different stacks. This argument works in 
the solid phase but not deep into the liquid phase where 
eventually stacks may lose their meaning, especially if 
the density of pancakes is large so the average distance is 
much smaller than A (high temperatures and high fields). 

Since the Josephson interaction is between nearest 
neighbor pancakes in adjacent planes it is quite straight 
forward to implement PBC. A pancake at the top plane 
(N p ) interacts with the closest pancake in the bottom 
plane 1 as well as with a pancake in plane N p — 1 . When 
calculating the lateral distance between pancakes we al- 
ways measure the "shortest distance" defined as follows: 
If the actual |Ax| separation is larger than L±/2 we sub- 
tract or add ~L\ depending on the sign of Ax, and sim- 
ilarly for Ay with L 2 replacing L\. This way the cor- 
rect distance is obtained even when the adjacent pan- 
cake in the plane above has exited the simulation cell 
and emerged close to the other side of the simulation 
cell. 

In the case of R 3> r g , string-string interactions 
that involve three and four-body interactions become 
important^. However near the melting transition for 
the range of magnetic fields investigated in this paper 
R « 0.25a ~ 1000A whereas r g « 5625A. Thus large 
transverse fluctuations for which the string-string inter- 
actions become important are statistically rare and can 
be neglected. Even in the case of tilted CDs we verified 
that the kinks present (see later discussion) are not fully 
developed Josephson strings. Thus in this work we do 
not include string-string interactions. 

The modeling of the interaction between pancakes 
and CDs was discussed by Goldschmidt and Cuansing 21 . 
There are two major sources of pinning: core pinning 
and EM pinning. Core pinning 2 arises when the vortex 
core overlaps with a normal state inclusion similar to the 
one inside a CD. Since condensation energy is lost in the 
vortex core, part or all of this energy is restored when a 
vortex resides inside a CD. EM pinning arise a 37 ' 38 when 
the supercurrent pattern around the vortex is disturbed 
by the non-conducting defect. These two mechanisms 
combine together to yield the expression for the poten- 
tial energy at a distance R away from the CD as felt by 
an individual pancake 2 -^: 
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where r r is the radius of the CD. The long range tail con- 
tribution rs —eodr^/R 2 is due mainly to the EM pinning 
and the short range flat region of depth « eod is due to 
the combined effects of EM and core pinning. This for- 
mula is valid provided r r > y2£. If r r < \f2^ a slightly 
different expression should be used, but we have deter- 



mined, see below, that in order to get good agreement 
with experiments we need to use r r > y/2£. Note that in 
this expression eo is temperature dependent and is pro- 
portional to 1 — T/T c . For example the pinning strength 
reduces to realistic values of 113 K at T=80 K and 79K 
at T=83K assuming T c = 90K. 

After experimenting with CDs of different radii (10 nm 
< r r < 30 nm) we concluded that in order to obtain a 
good agreement with experimental results we needed to 
choose r r w 30 nm in which case the formula given by 
Eq. (|2.6|) is the correct one. Different heavy ions with 
different energies give rise to tracks of various sizes, nor- 
mally reported to be in the range of 4 — 20 nm in di- 
ameter. It may be that the actual damage caused by a 
track exceeds its apparent size, thus leading to a higher 
effective radius. For tilted CDs it may be appropriate 
to take a potential with an elliptical cross section rather 
than a circular one. Since our simulations for tilted CDs 
are mainly for angles of 45° which are not that large and 
since the radius of our CDs is already larger than the ac- 
tual experimental radius, we did not use elliptical cross 
sections in the present simulations. 



C. Details of the Simulations 

The simulation cell is divided into a 800 x 692 mesh of 
small cells of area h x h each, where h = J Nfi / 800 (in 
units of ao). We tabulate the functions Go and Goc using 
a less refined cell structure creating two large 200 x 174 
matrices. During the simulations we use the tables as a 
lookup to calculate the pair interaction and use interpo- 
lation to find the applicable values. For each table we 
also calculate the negative of the gradient and save the 
two components of the gradient in their own tables. We 
also tabulate the Josephson interaction and its gradient 
and the columnar defect potential and its gradient but 
for these quantities we used a 800 x 692 matrices. When 
simulating we allow pancakes to move to arbitrary real 
locations, but in order to calculate the forces we divide 
the actual position by h and round to the nearest integer 
to use for the lookup tables. 

In each simulation step we move all the pancakes at 
the same time, using the instantaneous forces. This is 
done using a time step At. In these simulations we used 
the second order Runge Kutta method to advance the 
solution of the differential equations, which requires an 
additional virtual move for any real move. It is very im- 
portant to chose the time step correctly, especially when 
CDs are present. We used a variable size time step that 
limit how much a pancake can move in each step. The 
maximum distance allowed was 0.06 ao- Thus if the max- 
imum force on the pancakes is very large, a small time 
step is used. 

In the simulation we apply a procedure to implement 
flux cutting and recombination. We assume that within 
the coupled-planes model, vortices may switch connec- 
tions to lower their elastic energy (in this case Josephson 
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energy) when they cross each other—. In the simulation 
we construct two matrices of size Nfi x N p which we call 
the "up" matrix and the "down" matrix. For a given 
pancake i in plane p the "up" matrix points to the pan- 
cake in the plane p + 1 (or 1 if p = N p ) that is connected 
to the given pancake (i,p) via a Josephson interaction. 
Generally this is the pancake closest to the given pan- 
cake in the next plane. The "down" matrix similarly 
points to the closest pancake below (or in plane N p for 
p = 1). When we start from an initial configuration in 
which the FLs are straight stacks of pancakes, the ma- 
trices simply point to the pancake just above or below a 
given pancake. When constructing the force matrix after 
each time step we check if indeed the "up" matrix points 
to the closest pancake above. If there is a closer pan- 
cake than the one given by the pointer then we find out 
its parent in the plane p by using the down matrix, and 
we check if switching the two connections will decrease 
the sum of the squares of the two distances. If it does 
we cut and switch connections and update the "up" and 
"down" matrices. We term this process an "exchange". 
The reason we use the square of the distances is that in 
most instances the Josephson interaction is proportional 
to the square of the transverse distance. This procedure 
mimics the actual dynamics in which we expect the mag- 
netic flux to choose a path that minimizes the Josephson 
energy. Wc implement the flux cutting procedure after 
every update move of the system, but not during the 
virtual half-step in the Runge-Kutta procedure. 

Our simulations show that just below the melting tran- 
sition, even though some exchanges occur, they soon re- 
verse or undo themselves, and thus they do not lead to 
what we refer to as an entangled state where composite 
loops or permutations are abundant. On the other hand 
when exchanges proliferate through the system, a phe- 
nomenon that occurs in our simulations just above the 
melting transition, the exchanges do not undo each other, 
and the system of FLs changes from being composed of 
simple loops each made up of a single FL ending on it- 
self, to a system composed mainly of composite loops 
that wind up several times around the simulation cell 
in the z direction before returning to the original point. 
The reason that the exchanges proliferate above the melt- 
ing transition is that the transverse fluctuations become 
strong enough to overcome the potential barriers due to 
the repulsion among pancakes residing in the same plane 
and thus the crossing of FLs occur. 



D. Measured quantities 

We measured the following physical quantities. For 
details the reader is referred to our earlier wor k 18 ' 33 . The 
translational structure factor was also measured but we 
do not report on it in the present paper. 



1. Energy 

We calculated the average energy per pancake. The 
total energy is the sum of the contributions from the 
EM energy of all pairs of pancakes, the Josephson energy 
of nearest neighbor pancakes in adjacent planes and the 
binding energy of pancakes trapped by columnar defects. 
For the first part, the EM energy of a perfect Abrikosov 
lattice at the same temperature and magnetic field was 
subtracted from the total EM energy^. 

2. Mean square deviations 

For each individual flux-line we define the position of 
the lateral center of mass as Rcm = R(i,p)/N p , where 
the sum goes over all the pancakes belonging to it. We 
then define the mean square deviations as 

Rf = (y^X R (j,p) ~ Rcm) 2 ) , (2.7) 

Where the sum is over all pancakes belonging to an in- 
dividual flux-line and the average is over all flux-lines of 
the system and then taking a time average. The melt- 
ing transition is expected to occur when this quantity 
satisfies 

Rf/ao > c L , (2.8) 
where cl is the Lindemann coefficient. 

3. Line entanglement 

As we allow for flux cutting and recombination, we can 
define the number N e /Nfi as that fraction of the total 
number of FLs which belong to loops that are bigger 
than the size of a "simple" loop. A simple loop is defined 
as a set of N p beads connected end to end (due to the 
periodic boundary conditions in the z direction). Loops 
of size 2N P , 3N p ... start proliferating at and above the 
melting temperature. 

4- drift velocity under an applied force 

For the dynamical measurements, a uniform driving 
force of proper magnitude was applied along positive x 
direction. The average drift velocity of each pancake was 
measured. Both the distribution and the average value 
of these velocities are studied. They provide evidence on 
the derealization transition and the melting transition. 

5. Parameters 

Parameters for BSCCO were taken as follows: Ao = 
1700A, d = 15A, £ = 30A and T c = 90 K. The tempera- 
ture dependence of A in this work was taken to follow the 
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Ginzburg-Landau convention A 2 (T) = X^/(1 — T/T c ) and 
similarly for £. See discussion in Ref. [33| on the agree- 
ment of this choice with experiments. For the anisotropy 
7 we have used value of 375. For columnar defects we 
have used a radius of 30nm. 



III. VERTICAL COLUMNAR DEFECTS 

In this section, we concentrate on the simulations of 
systems with vertical CDs. Both static measurements 
and dynamical measurements have been done. For the 
static measurements, we start from a hexagonal lattice 
of straight and vertical FLs. The CDs attract those FLs 
which are near them and rapidly become occupied. For 
temperatures close to and above the melting transition, 
the FLs also become entangled. It takes relatively long 
for the entangling process to reach equilibrium, especially 
near the transition. We took 50 units of model time 
(time is measured in unit of 770^60) for equilibration and 
another 50 units of model time on measurements. Sim- 
ulations have been done on systems with 36 FLs and 36 
layers. Columnar defect density was set to be 50G. 
Magnetic fields from 50G to 200G were investigated. For 
each field value, 10 different CD configurations were sim- 
ulated and the measurement results were averaged. 

According to the periodic boundary condition along z 
direction, the top of each FL is connected via the Joseph- 
son interaction with the nearest pancake in the bottom 
plane. Below the melting transition, the pancake usu- 
ally belongs to the same FL as the top pancake. We say, 
the FL form a simple loop. As temperature increases, the 
thermal motion of the pancakes starts to make some pairs 
of FLs get close enough that flux cutting and reconnec- 
tion happens, and some FLs become part of a non-simple 
loop. Figure 



la 



shows the fraction of FLs belonging to 
non-simple loops. Below the melting transition, none of 
the FLs entangle with others and the curves lie exactly 
at zero; above the delocalization transition, nearly uni- 
form liquid is formed and the curves approach one. The 
criteria for melting can be taken to be that the fraction 
of non-simple loops reaches « 0.1 - 0.2. It shows that 
the melting tempe ratur e decreases as the magnetic field 
increases. Figure 1(b) shows the average mean square 



deviation of the pancakes on each FL relative to their 
"center of mass"—. The results are in units of a§ for B 
= 100G. The melting transition is characterized by the 
sharp change of slope. For example, for B = 50G, the 
change happens around effective temperature 600, which 
is about 78K, with mean square deviation « 0.1 - 0.2. 

Another quantity of interest is the number of trapped 
pancakes. Usually, each CD is occupied by a single FL. 
The trapped FLs form a cage for the interstitial FLs. Be- 
low the melting, the interstitial FLs form crystallites in 
the pores. At temperatures slightly above the melting 
transition, the crystallites melt locally and form nanoliq- 
uid in the pores. Higher temperature enables the trapped 
FLs to escape out of the CDs and the cage also melt, so 
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FIG. 1: (color online) Static simulations for systems of 36 
layers with = 50G and different values of B. The curves are 
average over 10 CD realizations, (a) Fraction of FLs belonging 
to non-simple loops, (b) Mean square deviations of FLs. 



a homogeneous liquid is formed. To compare results of 
different magnetic fields, Fig. HJa) plots the fraction of 
occupied CDs, which is given by the fraction of trapped 
pancakes multiplied by FL to CD ratio. As a common 
character for all the field values, the fractional values de- 
crease as temperature is increased. The curve for B = 
50G is significantly lower than that of higher field values 
at low temperatures. Even though the FL to CD ratio is 
one, the repulsion between the FLs tries to make them 
evenly spaced and some FLs have to be interstitial and 
some defects become unoccupied. Figure [2fb) shows the 
normalized binding energy, which is the binding energy 
per plane per CD. As a consequence of Fig. [Ua), the 
values increase as the temperature increases because the 
fraction of occupied CDs decrease. As a comparison with 
our previous work>2i, Fig. [DJc) is the average binding en- 
ergy per pancake. 

In order to verify that our system is large enough and 
finite size effects will not significantly change the results, 
we carried out simulations with 64 FLs in 36 planes. A 
comparison between the results of 36 and 64 FLs is de- 
picted in Fig. [3] We see that the results are very similar, 
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(b) 

FIG. 2: (color online) Static simulations for systems of 36 
layers with = 50G and different values of B. The curves 
are average over 10 CD realizations, (a) Fraction of occupied 
CDs. (b) Normalized binding energy, (c) Average binding 
energy per pancake, (b) and (c) used the colors as in (a) to 
represent different magnetic fields. 



the transition for 64 lines is slightly sharper but this is 
not a significant effect. 

Next, we discuss the dynamical measurements which 
were performed in order to locate the derealization tran- 
sition and the irreversibility line. The system is also 
started from a hexagonal lattice of straight and vertical 
FLs. After 20 to 50 units of model time for equilibration 
without driving force, a uniform in-plane driving force 
/l is applied on all the pancakes along the x direction. 
We used Jl = 0.025 and 0.1, which are small compared 
with the binding force on the trapped pancakes. Five 
to thirty units of model time were used on equilibration 
with driving force and 20 to 50 units of model time were 
used on measurements. 

Under a driving force, the interstitial pancakes are 
much easier to drift relative to the trapped ones. For 
driving force /z, ^0.1 and below the derealization tran- 
sition, most of the trapped pancakes remain trapped dur- 
ing the time span of our simulation. The average velocity 
of each pancake was measured. Figure 0] shows its dis- 
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T/(1 -T/T )(K) 
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FIG. 3: (color online) Static simulations for systems with 36 
FLs (blue curves with circles) and 64 FLs (red curves). The 
curves are average over 10 CD realizations. B = 100G and 
Bj, — 50G. (a) Fraction of non-simple loop FLs. (b) Mean 
square deviations of FLs. 
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FIG. 4: Normalized distribution of average velocity along x 
direction. B = 100G, B^ — 50G, simulation of system with 
36 FLs and 200 layers. The driving force ft, was taken to be 
0.1 and 4 different CD configurations were done, (a), (b) and 
(c) are the distributions at temperatures in the vicinity of the 
derealization transition. 
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FIG. 5: Probability density of zero average velocity at vari- 
ous temperatures. Results for three different fields are shown 
here. Simulations on systems with 36 FLs and 36 layers. 
= 50G, driving force /l = 0.1. Average over 10 different CD 
configurations. 
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FIG. 6: Average drift velocity over driving force. Results 
with magnetic fields of 50G, 100G and 200G are shown. Used 
the same data as in Fig. [5] for B = 50 G and B = 100 G, 
where, f L = 0.1. For B = 200 G, f L = 0.025 was used. 



tribution for temperatures near the delocalization tran- 
sition. The distribution is normalized so that the area 
of the region under the curve is one. Below the transi- 
tion temperature, the distribution is characterized by a 
sharp peak at zero and a wider and lower peak centered 
above zero. The sharp peak corresponds to the trapped 
pancakes and the wide peak corresponds to the intersti- 
tial ones. The driving force we used was small enough 
that the trapped pancakes would not be pulled out of the 
columnar defects below the delocalization transition but 
large enough that the interstitial pancakes could move 
through the skeleton formed by trapped FLs. As tem- 
perature increases, the height of the sharp peak decreases 
and the wide peak becomes wider. The sharp peak dis- 
appears above the delocalization temperature and the 
center of the remaining Gaussian shaped wide peak ap- 
proaches 0.1, which is the value of the driving force. 

One should notice that the width and height of the 
distribution in Fig. 2] depends on the length of the mea- 
surement. The longer the measurement is, the narrower 
and taller the wide peak is. For random walk, the average 
distance the walkers travel is proportional to the square 
root of time, thus the average velocity of the walkers de- 
creases as the time of averaging becomes longer. The 
movement of the interstitial pancakes are different from 
a simple random walk. They are confined by other FLs, 
especially by the trapped ones. They are also confined 
by pancakes in adjacent layers through Josephson inter- 
action. However, our simulations do show that the width 
of the wide peak times the square root of the time span 
of the measurement is roughly constant near and above 
the delocalization transition. 

In Fig. the probability density of zero average ve- 
locity, i.e. the value of the average velocity distribution 
at zero, is plotted. Results for magnetic fields of 50G, 
100G and 200G are shown. The delocalization transition 
is characterized by the abrupt change of the slop of the 



plotted curves, which happens around effective temper- 
atures 600, 500 and 400 respectively, corresponding to 
78K, 76K and 73K. The delocalization transition tem- 
perature decreases as the magnetic field increases. 

It is possible to locate the irreversibility line from the 
time averaged drift velocities, by taking an average over 
all pancakes. The average velocity of the pancakes is de- 
picted in Figj6] For temperatures below the irreversibil- 
ity line, which is the threshold of mobility as a function 
of temperature, the curves are flat and the drift veloc- 
ity is small. Then at a temperature about equal to the 
melting temperature, their slopes increase sharply. The 
sharp rise in mobility is most clear for the B — 50G case, 
for which the melting transition and the delocalization 
transition happen at the same temperature. At low tem- 
peratures, the FLs are pinned collectively by the CDs. 
They are nearly immobile under a small driving force. 
When the temperature is raised high enough, that the 
interstitial pancakes can move easily in the pores, a small 
driving force will make the FLs drift, signaling the on- 
set of depinning. Above the delocalization temperature, 
the trapped pancakes also start to drift under the small 
driving force, making the ratio of average drift velocity 
over driving force approach one. 

The results above are summarized into a phase dia- 
gram (Fig[7]). It is divided into three regions by the 
melting line and the delocalization line Bdi- Be- 
tween these two lines, the vortex system is in a nanoliq- 
uid phase, in which the interstitial FLs are entangled and 
melted but are separated from the trapped FLs. To the 
left of the dashed line associated with the melting line, 
the system is in the nanocrystaline phase, in which all the 
FLs are separated from each other, the interstitial ones 
form triangular lattices in the pores surrounded by the 
trapped ones. The error bars indicate the width of the 
transition region due to the finite size of our system. The 
melting line B^ of a pristine system is also plotted for 
comparison. To the right of the dashed line associated 
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FIG. 7: (color online) Phase diagram in T-B plane. The melt- 
ing transition B^ of pristine system, the melting transition 
of system with columnar defects [B^ = 50 G) and the 
derealization transition Bdi of system with columnar defects 
(B^ = 50 G) are shown. 



with the delocalization line, the system is in a homoge- 
neous liquid phase, where the effect of the CDs becomes 
less important and all the FLs become entangled. 



IV. TILTED COLUMNAR DEFECTS 

It is interesting to study the effect of the tilting of 
CDs on the properties of the vortex system. Hwa, Nel- 
son and Vinokur— , Nelson and Vinokur 22 ^, and more 
recently Refael et alM studied the effect of a tilted mag- 
netic field on the properties of the Bose glass system. In 
their investigation the CDs are along the c-axis and the 
magnetic field has a component H±_ along the a-b planes, 
thus forming an angle 8 between the direction of the ap- 
plied field and the CDs. When the angle 8 is small the 
vortices which are localized on the CDs are completely 
trapped along the CDs and the Bose glass phase is pre- 
served. At larger 9 > 8 C (H± > Hj_), the vortices form a 
staircase structure hopping from one CD to the next. For 
different vortices those kinks tend to align in "chains". 
For an even larger angle 8 > 8 a the vortices follow the 
field direction and are essentially unaffected by the corre- 
lated nature of the CDs. This picture appears to describe 
correctly the situation in YBCO in the Bose glass regime 
(B < B,),) and is supported by various experiments, both 
for thermodynamical 4 ^ and dynamical propertie a 45 ' 46 . 

For BSCCO which is much more anisotropic the situ- 
ation appears more complicated and less clea r 47 ' 48 . Fur- 
thermore in this case there are interesting results on the 
nanosolid phase which exists for B > B^ and is differ- 
ent from the Bose glass phase. There are no theoretical 
investigations on how the nanosolid phase and its two 
stage melting transition might be affected by a tilting 
angle between the direction of the CDs and the magnetic 
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FIG. 8: (color online) Static measurements on systems with 
vertical CDs (blue curves with circles) and titled CDs (red 
curves, CDs are tilted by 45° away from c axis), (a) Fraction 
of occupied CDs. (b) Total energy per pancake, (c) Mean 
square deviation of FLs. (d) Exchange rate, (e) Fraction of 
FLs belonging to non-simple loops. 



field. To try to answer this question we recently car- 
ried out experiments and numerical simulations to try to 
understand the properties of such a system 2 -. Here we 
give more details on the method and the results of the 
simulations. 

In this work we limited our investigation to the case of 
a vertical external magnetic field (parallel to the c axis), 
and considered only the case of tilted CDs at various 
angles. Simulations with 36 layers are not enough any 
more, since the shift of the top of the CDs relative to their 
bottom would be too small. We used 200 layers, which is 
still manageable by the current computing power. Most 
of the simulations were done on systems with 36 FLs 
and 18 CDs. The magnetic field was set to be 100G. For 
comparison, both systems with vertical CDs and systems 
with tilted CDs were simulated. The comparison between 
the two systems was done for the same value of B^ — 
B§ cos 8, in order to match the experimental situation 
in which the ratio of vortices to CDs remains fixed as 
8 changes 2 ^. Both static and dynamical measurements 
were carried out and the results are averaged over four 
different random CD configurations. 

Figure[S]shows the results of static measurements. Fig- 
ure [5£a) shows the fraction of occupied CDs, i.e. the 
number of trapped pancakes over the number of trap- 
ping sites. Over the temperature range of our simula- 
tions, the vertical case and the tilted case give the same 
results within the error bars. At effective temperatures 
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below 350, the ratio is close to one, so the CDs are al- 
most fully occupied. In Fig. |5fb), the total energy per 
pancake are almost the same in the two cases. In Fig. 
HJc), the mean square deviation of the FLs in the tilted 
case is a little higher than that of the vertical case at 
low temperatures. Since our measurements are relative 
to the center of mass of each FL, and as will be shown 
later, the FLs are tilted by the tilted CDs, the deviations 
are raised to higher values. The sudden change in the 
curves' slopes reflects the position of the melting transi- 
tion of the interstitial FLs. Figure [5Jc) might suggest a 
one degree increase in the melting temperature for tilted 
CDs, but this is within the simulation error. Figure ^d) 
shows the number of "exchanges" per unit time, i.e. the 
rate of cutting and recombining of FLs. This number 
is under 1000 below the melting transition and increases 
rapidly above the melting. We see that it is the same 
for both straight and tilted CDs. Figure [5Je) shows the 
amount of entanglement, i.e. the fraction of non-simple 
loops. The flat area of the curves is an indication for the 
melting transition. First, we note that for straight CDs, 
the fraction of non-simple loops at the transition for 200 
layers is higher than the corresponding value for 36 lay- 
ers, which is about 0.2. This is because for longer FLs, 
there is more chance to cross and also to form occasional 
kinks involving the CDs. This figure is the only one that 
shows a significant difference between the straight and 
tilted CDs in terms of the absolute value of their entan- 
glement. This is likely due to the fact that FLs trapped 
on CDs are tilted and so are those FLs in the vicinity of 
a CD. Thus the top and bottom of those FLs are not ad- 
jacent and they may connect to other FLs to form loops. 
This interpretation is favored by the fact that the num- 
ber of exchanges for the tilted CD case is not larger than 
for the case of straight CDs but the amount of their en- 
tanglement is larger nevertheless. Thus this criterion (of 
entanglement) using the boundary conditions may not be 
a good indication for the melting transition in the case of 
tilted CDs as it is for the case of pure system or a system 
with straight CDs. The flat plateau region of the curve 
should give a rough indication where the melting takes 
place. Thus for tilted CDs, we base our criterion for the 
melting on the transverse fluctuations of the FLs and on 
the number of exchanges per unit time rather than on 
the amount of entanglement. 

We now account for the energies involved: The elec- 
tromagnetic energy cost per layer associated with the 
tilting of an infinitely long stack of pancakes is given 
byi^ eodln((l + cos0)/2cos0). This cost is exceeded, 
for angles less than about 75°, by the energy gain from 
the pinning of columnar defects, which is of the order of 
eod per layer. In addition, one should also include the 
energy B • H/(47r), which favors to align the average 
vortex orientation with H , which is implemented in the 
simulation by the periodic boundary conditions: When 
a stack is tilted, the stacks of images along the z direc- 
tion are not positioned along the tilted direction but are 
above (and below) the center of mass of the tilted stack 



and because of the electromagnetic interaction they pull 
the stack in the z-direction, acting like the B • H/(47r) 
term. The competition between the binding energy gain 
and the loss of electromagnetic energy associated with a 
tilt favors the formation of kinks. The trapped FLs tilt to 
take advantage of the binding energy of the tilted CDs, 
and the interstitial FLs tilt because of repulsion by the 
trapped pancakes. But in order to relieve the frustra- 
tion of the electromagnetic interaction, kinks are formed 
and the FLs jump back to align with the z axis which is 
the external field direction. The larger the tilting angle 
of the CDs the more kinks appear along the FLs. The 
kinks on neighboring FLs tend to align with each other. 
The formation of the kinks increases the Josephson in- 
teraction part, however, in high anisotropic materials like 
BSCCO, this cost is negligible compared with the total 
energy 49 . Near the transition the transverse fluctuations 
due to the tilt-kink combination is not much different 
from the transverse fluctuations due to temperature fluc- 
tuations. Our simulations show that the vortices form a 
tilted rigid matrix which is parallel to the CDs, and the 
CDs are almost fully decorated by pancakes. Because of 
the kinks, the interstitial vortices can switch cage after 
a certain distance along the z direction, but each section 
is still caged by the trapped pancakes on the CDs. Also 
interstitial FLs can themselves by trapped on CDs for 
some portion of their length. In the vicinity of the melt- 
ing transition (T — 73K to 77K), considering the error 
bars of the measurements, we found that systems with 
CDs tilted at 45° have nearly the same total energy as 
those with non-tilted CDs. One can conclude that the en- 
hanced caging potential created by the CDs, which raises 
the melting temperature^, is not significantly affected 
when increasing the angle between B and the CDs, thus 
the thermodynamic properties, i.e., the melting transi- 
tion and the delocalization, are angle independent. This 
has been shown in our static simulations in Fig. [5J and 
will be further demonstrated for the delocalization tran- 
sition in the following discussion. 

The effect of tilting of the CDs on the delocalization 
transition is also investigated. Figure [9] shows the prob- 
ability density of zero average velocity for both vertical 
case and titled case. The two curves coincide within the 
range of the error bars, especially at temperatures near 
the delocalization transition. The average velocity dis- 
tributions of these two cases are also compared and one 
cannot see any significant difference. One can conclude 
that the tilting of CDs does not have noticeable effect on 
the delocalization transition of the vortex system. 

Snapshots are shown here on systems without a driving 
force. FigurefTOTa') andfTOTb - ) shows the top view of a sys- 
tem with CDs tilted by 45° at 60K and 75K respectively. 
At 60K, the system is in the nanosolid phase and indi- 
vidual FLs can be distinguished clearly. The interstitial 
FLs try to form hexagonal lattices. At 75K, the system 
is in nanoliquid phase. Interstitial FLs become entangled 
and the pancake vortices distribute nearly uniformly in 
the pores. But the temperature is not high enough to 
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FIG. 9: (color online) Probability density of zero average 
velocity, the blue curve with circles is for vertical CDs and 
the red curve with triangles is for tilted CDs. Driving force 
h = 0.1. 



(a) T = 60K 



(b) T = 75K 
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FIG. 10: (color online) Snapshots of vortex stacks and CDs 
tilted at 45° (a, b and c) and 80°, Bf = 36 G. (d) for B = 100 
G. Pancakes belonging to the same stack are connected. Free 
pancakes are in blue and trapped ones in red. Only CDs at 
the bottom layer are shown (green). (a)(b) Projection onto 
a-b plane (top- view), (a) Nanosolid phase, (b) Nanoliquid 
phase, (c) Projection of (a) onto a-c plane (side view), first 
row is shown, (d) Side view with CDs at 80°. 



FIG. 11: (color online) Snapshot of a 2D system with 200 
planes, 6 FLs and 3 CDs. Pancakes belonging to the same 
stack are connected. Trapped ones are in red and interstitial 
ones are in blue. CDs are evenly spaced and titled at 73° so 
that they are continuous at the boundaries. T = 70 K. 



into segments, which are titled to some degree along the 
direction of the CDs. Because of the kinks, the FL seg- 
ments jump back periodically, keeping the FLs vertical on 
average. At higher temperatures, The kinked structures 
and tilting of interstitial FLs still exist but are compli- 
cated by the entanglement of the FLs. In Fig. [TOM), 
the CDs are tilted by a larger angle than in Fig. fTUTc). 
as a consequence, more kinks appear and the separation 
between the trains of kinks becomes smaller. 

Under large enough driving force and temperature in 
the proper range, FL kinks can move along z direc- 
tion and the pancake vortices drift along the direction 
of the driving force. To show this, we simulated a two- 
dimensional system with six FLs and three titled CDs 
(Fig. Hip . The pancakes were confined to move only in 
the x-z plane. The CDs were arranged such that they 
are continues at the boundaries. We initialized the FLs 
to be vertical, straight and evenly spaced. After simu- 
lating for enough time to the let the system equilibrate 
without a driving force, two trains of kinks were formed. 
A driving force of magnitude fa = 0.5 was then applied 
along the positive x direction. At 70 K, the kinks moved 
downward altogether with a constant speed of about 13 
Cu02 planes for every ten time units. 



make trapped pancakes delocalize yet. The nanoliquid 
is repelled by the trapped FLs, thus forming void spaces 
surrounding the CDs. Figure [TUT c) is a side view of the 
first row of the FLs in Fig. [TUT a). Usually, a CD is occu- 
pied by more than one FL and kinked structures emerge. 
The kinks go along the a-b plane and extend for less 
than ten CuC"2 layers. There are also kinks on intersti- 
tial FLs. They tend to reside at the same layers with 
those attached to CDs and form trains of kinks^i , which 
resemble the Josephson and Abrikosov crossing lattices 50 
found in BSCCO with field tilted with respect to the c 
axis. Another feature is that the interstitial FLs break 



V. CONCLUSIONS 

To conclude we have used Langevin molecular dynam- 
ics simulations to investigate the melting transition of the 
vortex matter in the presence of vertical CDs and in the 
presence of tilted CDs. For vertical CDs our results are 
in full agreement with o ur p revious multilevel MC sim- 
ulations reported in Ref.[21]. Furthermore we measured 
also dynamical properties of the system under an exter- 
nal current that applies force on the vortices and these 
located the melting and derealization transitions at the 
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same temperatures as the static measurements, giving a 
more precise identification of the location of the delocal- 
ization transition. For the melting transition the dynami- 
cal measurement actually identifies the irreversibility line 
associated with partial depinning, and this occurs at the 
same temperature as the melting for the case of vertical 
CDs. 

For the case of tilted CDs, the porous vortex matter 
in BSCCO is shown to preserve its thermodynamic prop- 
erties even when the field is tilted away from the CDs. 
The positions of the melting and derealization lines re- 
main unchanged while the tilting angle of the CDs vary, 
at least up to angles of 45°. Our simulations show that 
increasing the angle between the field and the CDs leads 
to formation of weakly pinned vortex kinks, while pre- 
serving the basic structure of a rigid matrix of pancakes 
residing along the CDs with nanocrystals of interstitial 
vortices embedded within the pores of the matrix. We 
also showed that under the influence of an applied force 
(like a current) the kinks slide down and lead to a global 
motion of the pancakes in the direction of the applied 



force. Thus we expect that for large tilting angle between 
the CDs and the magnetic field, which result in the pro- 
liferation of kinks, the irreversibility line should move to 
lower temperature, closer to the location of the melting in 
the pure system (no CDs). To measure this shift explic- 
itly in the simulations as a function of the angle requires 
much more computational time than currently available, 
but it was verified in the recent experiments^. 
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